This set of simulations aims to compare different workflows for arrranging healthcare workers at a blood donation site to minimize risk of viral traqnsmission. We are using it to suggest workflow improvements for blood donation sites in the context of COVID-19. This model DOES NOT aim to accurately estimate the risk of healthcare workers or blood donors contracting COVID-19. It is just a tool used to compare workflow designs.
We are going to work with a very simpified model of virus transmission where the risk of virus spread from any given interaction is a combination of 3 factors:
proportion_positive <- 0.1
transmission_rate <- 0.2
ppe_failure_rate <- 0.1
interaction_risk <- proportion_positive * transmission_rate * ppe_failure_rate
cat(interaction_risk,"or",interaction_risk*100,"%")
## 0.002 or 0.2 %
Therefore, all three factors can be combined into the single relevant parameter interaction_risk. We will run simulations across a range of the parameter interaction_risk so tat we can encompass the uncertainty in all three of the factors which define it.
The basic model will be this: a donor enters and must complete a 1) health screening and 2) blood collection. Each of those actions involves an interaction with a healthcare worker.
In each case, there is a set number of healthcare workers who can either act as screeners (S) who do the healthcare screening or phlebotomists (P) who draw the blood. The donor sees one screener and one phlebotomist.
The simulation will start with all healthcare workers being COVID-negative and simulate their random interactions with donors for the specified number of days. The output will be how many of the healthcare workers have contracted COVID-19 after X amount of days under various conditions for interaction_risk.
We will start with the following two necessary parameters:
We are going to start with our simplest workflow design: the blood donation center has 10 healthcare workers, all of whom can be S or P. They interact with donors randomly. The simulation is carried out by the function simulate_hcw_risk.
hcw_risk_vec <- simulate_hcw_risk(n_hcw = 6,
n_donors = 20,
n_days = 30,
interaction_risk = 0.002)
We can plot the proportion of healthcare workers contracting the virus over time.
Since we are doing a random simulation for each interaction, we should simulate 100 times and take the average.
n_sim <- 100
hcw_risk_mat <- matrix(nrow = n_sim, ncol = 30)
for (i in 1:n_sim){
hcw_risk_mat[i,] <- simulate_hcw_risk(n_hcw = 6,
n_donors = 20,
n_days = 30,
interaction_risk = 0.002)
}
hcw_risk_mean <- apply(hcw_risk_mat, MARGIN = 2, FUN = mean)
And we can plot the mean value over the 30 days.
We are also interested in how the transmission is affected by the interaction risk, which is a lumped parameter which describes how likely a healthcare worker is to contract COVID from interacting with a blood donor. Below, we will vary the interaction risk, simulate for a given number of days 100 times, and take the average at the final timepoint as the final proportion of healthcare workers contracting the disease.
# Set which values of interaction_risk to explore
interaction_risk_vec <- seq(from = 0.0001, to = 0.02, by = 0.001)
# Set the number of simulations
n_sim <- 100
# Initialize a vector
interaction_risk_results <- rep(0, times = length(interaction_risk_vec))
# Simulate for each value of interaction_risk
for (i in 1:length(interaction_risk_vec)){
# Initialize matrix for results
hcw_risk_mat <- matrix(nrow = n_sim, ncol = 30)
# Simulate the model across chosen number of simulations
for (j in 1:n_sim){
hcw_risk_mat[j,] <- simulate_hcw_risk(n_hcw = 6,
n_donors = 20,
n_days = 30,
interaction_risk = interaction_risk_vec[i])
}
# Take the average across all simulations
hcw_risk_mean <- apply(hcw_risk_mat, MARGIN = 2, FUN = mean)
# Store into result vector
interaction_risk_results[i] <- hcw_risk_mean[length(hcw_risk_mean)]
}
And we can plot the result:
Within the same basic model, we want to include the possibility that once infected, a healthcare worker can become an asymptomatic carrier and spread the disease to donors they interact with in the future.
This can be achieved by using the same simulation but incorporating a risk of interaction from healthcare worker to donor. The risk of this happening (donor_interaction_risk) can be calculated similarly to the interaction_risk above. There is no proportion_positive term since the chance of this spread only happens if the healthcare worker is positive for the virus.
transmission_rate <- 0.2
ppe_failure_rate <- 0.1
donor_interaction_risk <- transmission_rate * ppe_failure_rate
cat(donor_interaction_risk,"or",donor_interaction_risk*100,"%")
## 0.02 or 2 %
The simulation takes place within the function simulate_donor_risk:
donor_new_cases <- simulate_donor_risk(n_hcw = 6,
n_donors = 20,
n_days = 30,
interaction_risk = 0.002,
proportion_positive = 0.1)
We can plot the result of new cases over the timecourse:
Again, we can simulate 100 times and take the mean, then plot over 30 days.
n_sim <- 100
donor_risk_mat <- matrix(nrow = n_sim, ncol = 30)
for (i in 1:n_sim){
donor_risk_mat[i,] <- simulate_donor_risk(n_hcw = 6,
n_donors = 20,
n_days = 30,
interaction_risk = 0.002,
proportion_positive = 0.1)
}
donor_risk_mean <- apply(donor_risk_mat, MARGIN = 2, FUN = mean)
donor_risk_mean_df <- data.frame(donor_risk_mean)
ggplot(data = donor_risk_mean_df, aes(y = donor_risk_mean, x = 1:length(donor_risk_mean))) +
labs(x = "Days", y = "", title = "Number of Donors Contracting Virus from Donation Site",
subtitle = paste("Average across",n_sim,"simulations")) +
geom_line() + theme_minimal() +
theme(plot.title = element_text(size = 18), axis.title.y = element_text(size = 16))
And we can also explore the impact of the interaction_risk parameter on the number of new cases acquired from the blood donation site at the end of the time period.
We will compare the following 5 ways of aranging health screeners (S) and phlebotomists (P).
The simulations above represent case 3. Below we will repeat the same analysis steps for cases 1, 2, 4 and 5.
The proportion of healthcare workers infected over 30 days.
The number of new donor cases acquired at the donation center.
The proportion of healthcare workers infected over 30 days.
The number of new donor cases acquired at the donation center.
The proportion of healthcare workers infected over 30 days.
The number of new donor cases acquired at the donation center.
The proportion of healthcare workers infected over 30 days.
The number of new donor cases acquired at the donation center.
The proportion of healthcare workers infected over 30 days.
The number of new donor cases acquired at the donation center.
Comparing the 5 workflow strategy in the transmission to healthcare workers and to donors.
We can compare the proportion of healthcare workers who contract the virus over time. And how the final number of infected healthcare workers varies with interaction risk.
hcw_risk_mean_combined <- data.frame(case1 = hcw_risk_mean_case1,
case2 = hcw_risk_mean_case2,
case3 = hcw_risk_mean_case3,
case4 = hcw_risk_mean_case4,
case5 = hcw_risk_mean_case5,
time = 1:n_days)
hcw_risk_mean_combined_melted <- melt(hcw_risk_mean_combined, id = 'time')
ggplot(data = hcw_risk_mean_combined_melted, aes(y = value, x = time, color = variable)) +
labs(x = "Days", y = "", title = "Proportion of Healthcare Workers Contracting Virus",
subtitle = paste("Average across",n_sim,"simulations")) + ylim(0,1) +
geom_line(size = 1) + theme_minimal() +
theme(plot.title = element_text(size = 16), axis.title.y = element_text(size = 16))
# Varying interaction risk
interaction_risk_combined <- data.frame(case1 = interaction_risk_results_case1,
case2 = interaction_risk_results_case2,
case3 = interaction_risk_results_case3,
case4 = interaction_risk_results_case4,
case5 = interaction_risk_results_case5,
interaction = interaction_risk_vec)
interaction_risk_combined_melted <- melt(interaction_risk_combined, id = 'interaction')
ggplot(data = interaction_risk_combined_melted, aes(y = value, x = interaction, color = variable)) +
labs(x = "Interaction Risk", y = "", title ="Proportion of Healthcare Workers Infected by Final Day",
subtitle = paste("Averaged across",n_sim,"simulations")) + ylim(0,1) +
geom_line(size = 1) + geom_point() + theme_minimal() +
theme(plot.title = element_text(size = 16), axis.title.y = element_text(size = 16))
We can compare the average number of donors acquiring virus from the donation center over time. And then look at the total number of donation center acquired cases at the end of the timecourse, under different interaction risks.
donor_risk_mean_combined <- data.frame(case1 = donor_risk_mean_case1,
case2 = donor_risk_mean_case2,
case3 = donor_risk_mean_case3,
case4 = donor_risk_mean_case4,
case5 = donor_risk_mean_case5,
time = 1:n_days)
donor_risk_mean_combined_melted <- melt(donor_risk_mean_combined, id = 'time')
ggplot(data = donor_risk_mean_combined_melted, aes(y = value, x = time, color = variable)) +
labs(x = "Days", y = "", title = "Number of Donors Contracting Virus from Donation Site",
subtitle = paste("Average across",n_sim,"simulations")) +
geom_line(size = 1) + theme_minimal() +
theme(plot.title = element_text(size = 16), axis.title.y = element_text(size = 16))
# Varying interaction risk
interaction_risk_donor_combined <- data.frame(case1 = interaction_risk_results_donor_case1,
case2 = interaction_risk_results_donor_case2,
case3 = interaction_risk_results_donor_case3,
case4 = interaction_risk_results_donor_case4,
case5 = interaction_risk_results_donor_case5,
interaction = interaction_risk_vec)
interaction_risk_donor_combined_melted <- melt(interaction_risk_donor_combined, id = 'interaction')
ggplot(data = interaction_risk_donor_combined_melted, aes(y = value, x = interaction, color = variable)) +
labs(x = "Interaction Risk", y = "", title ="Number of New Cases Acquired from Donation Site at Final Day",
subtitle = paste("Averaged across",n_sim,"simulations")) +
geom_line(size = 1) + geom_point() + theme_minimal() +
theme(plot.title = element_text(size = 16), axis.title.y = element_text(size = 16))